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ABSTRACT 

We perform high resolution A^-body+SPH simulations of isolated Milky- Way-like 
galaxies and major mergers between them, to investigate the effect of feedback from 
both an active galactic nucleus (AGN) and supernovae on the galaxy's evolution. Sev- 
eral AGN methods from the literature are used independently and in conjunction with 
supernova feedback to isolate the most important factors of these feedback processes. 
We find that in isolated galaxies, supernovae dominate the suppression of star forma- 
tion but the star formation rate is unaffected by the presence of an AGN. In mergers 
the converse is true when models with strong AGN feedback are considered, shutting 
off star formation before a starburst can occur. AGN and supernovae simulated to- 
gether suppress star formation only slightly more than if they acted independently. 
This low-level interaction between the feedback processes is due to AGN feedback 
maintaining the temperature of a hot halo of gas formed by supernovae. For each 
of the feedback processes the heating temperature is the dominant parameter rather 
than the overall energy budget or timing of heating events. Finally, we find that the 
black hole mass is highly resolution dependent, with more massive black holes found 
in lower resolution simulations. 

Key words: galaxies: formation - galaxies: evolution - galaxies: active - black hole 
physics - methods: numerical - galaxies: interaction 



1 INTRODUCTION 

In recent years, it has been recognised that the feedback en- 
ergy associated with black hole accretion has an important 
influence on structure formation, on both galactic and ex- 
tragalactic scales (e.g. iMcNamara &: Nulsenll2007^ . Active 
gala ctic nuclei (AGN) are now ob served as far back as z ~ 7 
fe.g. iFan et"alll200ll : JKurk et al.. 2007; Mortlock 2012) and 
it is now thought that most, if not all galaxies host a super- 
massive black hole at their centre. Moreover, in systems 
where black hole mass can be estimated, it has been shown 
to correla te with the ho st galaxy's spheroid velocity disper- 
sion (e.g. iFerrarese fcM erritt 2000; Trcmainc ct al. 2002), 
luminosity and ste l lar mass (e.g. iMagorria n et al.l 1 19981: 
Haring fc RbJ 12004 iBennert et all [2011. ; .McConnefl fc Mai 
2013). These results point to a fundamental link between 
the formation and evolution of galaxies and their central 
blac k holes; a position supported by theoretical argument s 
fe.g. lSilk fc ReeJl998l : [King|2003l ; lishibashi fc Fabianl2012h . 
Further evidence for the importance of AGN in the evolu- 
tion of galaxies is provided by semi-analytic modelling of 
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galaxy formation ( e.g. Benson et al.ll2003l : lBower et al.ll2006l : 
ICroton et aflbOoAlGuo et al.ll201ll ) where AGN feedback is 
invoked to reproduce the sharp drop-olf in number density 
at the bright end of the galaxy luminosity function (see also 
lBinnevll200J '). 

The modell ing of AGN i n galaxy fo rmation simula- 
tions s et out by ISpringel et al.l ([2005) and Di Matteo et al.l 
(|2005l . I2OO8I ). showed for the first time that black hole 
self regulation can be achieved and the observed Mbh — 
AfbuigG relation reproduced within a fully cosmological sim- 
ulation. Thi s method i s now widely employed in sinm- 
lations (e.g. iRobertsoiTe t al. 2006; Di Matteo ct al. 200i; 
Khalatvan et al.ll2008l; I Johansson et al.ll2009l : I Fabian et al.l 



20ld : lDeGraf et al.ll2012l ). although is by no means the only 



implementation currently in use and many aspects of the 
model are questioned. One notable assumption is that the 
accretion o n to the black hole is giv e n by the Bondi-Hoyle- 
Lyttle t on jHovle fc LyttletonI Il939l : iBondi fc HwJi Il944l : 
iBondil 1 19521 ; 'Bondi') formula which assumes spherically 
symmetric accretion from gas which is stationary at infin- 
ity. This has been challenged by a number of authors, for 
example B ooth fc Schaye (2009 ) dispute the correction re- 
quired for under-resolved accretion flow and central densi- 
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ties around the black hole. Other authors further dispute 
the Bondi approach due to the inherent assumptions of 
spherical syrnmetry and neRliRible ang u lar momentu r n (e.g . 
iDebuhr et allboiol : iPower et al.]|201ll ). iHobbs et all tom 
argue that a proper accretion rate estimate should take into 
account the gravitational potential of all matter, not just 
the black hole. Each of these modifications, as well as oth- 
ers taking into account more com plex physics such as black 
hole spin (IFanidakis et al.ll201ll ). will yield different accre- 
tion rates onto the black hole and thus alter its final mass. 
Furthermore, the method by which AGN feedback energy is 
deposited in the surrounding gas is often varied with some 
work s arguing for thermal feedback (e.g. Di Matteo et al.l 
I2OO5I : iBooth fc Schav3 l2009l : ISiiacki fc Springell l200d) and 
others for kinetic or momentum feedback (e .g. iDebuhr et al.l 
l2010l : |Power et al.ll201ll : [Dubois et al.ll2012l ) on various phys- 
ical and numerical grounds. A key criticism of thermal feed- 
back is that much of the energy used in raising the gas tem- 
perature is lost quickly through radiative cooling, a problem 
overcome in the method of Booth & Schayc (2009) which as- 
sumes that energy is stored until a gas particle can be heated 
to a high temperature. 

Another issue that is now receiving attention is the 
interaction between feedback from AGN and supernovae. 
Typically, these two feedback mechanis ms are treated sep- 
arately in semi-analytic models ( e.g. iBower et al.l 120061 : 
ICroton et al. 20061: IGuo et al.ll201ll ). However, recent work 
bv lBooth fc Schavel (|2012l ) using cosmological hydrodynam- 
ical simulations, casts doubt on the validity of this assump- 
tion. Their findings indicate that the combination of the two 
feedback processes has a weaker effect on the star formation 
rate than expected if they were acting independently. Due to 
the uncertainty in our present modelling of AGN feedback, 
it may be premature to assume that this interplay between 
AGN and supernovae is a physical effect and is not merely 
an outcome of the specific combination of models employed 
bv lBooth fc Schavel (|2012l ). 

In this paper, we seek to elucidate this interplay by 
studying how AGN and supernova feedback (both inde- 
pendently and combined) alter the star formation rate 
within simulations of isolated and merging disc galax- 
ies. Performing such idealised simulations is a commonly 
employed method when studying such numerically de- 
manding processes, allowing results to be produced with 
higher resolution than is normally possible in full cos- 
mological simulations. Such systems have been used in 
the past to simula te a range of topics in cluding galac- 
tic dyna mics (e.g. Springel fc White! 1 1999 ). star forma- 
tion (e.g.lSchave fc Dalla Veccliiall2008l'), s upernova feedback 



(e.g. iDaUa Vecchia fc Schave f'2008'. '2 0121) and AGN imple - 



mentations (e.g. lDi Matteo e t al. 20051: iDebuhr et aLllJoTol ) 



Our main results focus on three methods used in the litera- 
ture for implementing black hole growth and AGN feedback, 
including a novel adaptation of the method suggested by 
I Power et al.l (|2011^ . To simplify the interpretation of our re- 
sults we adopt the same models for all other sub-gnd physics 
included within our simulations (including supernova feed- 
back, star formation and interstellar gas). 

The remainder of the paper is organised as follows. In 
Section [5] we describe our numerical method, namely how 
we set up our initial conditions and give details of the vari- 
ous sub-grid implementations, including the supernova and 



AGN feedback models. Our main results are presented for 
isolated galaxies in Section |3] and merging galaxies in Sec- 
tion 131 while in Section [5] we take a closer look at some of 
our model assumptions. Finally, conclusions are drawn and 
future work outlined, in Section (6] 



2 NUMERICAL METHOD 

All simulations discussed in this paper are performed us- 
ing a modified version of the publicly-available A^-body -I- 
Smoothed Particle H ydrodynamics (SP H; see e.g.lMonaghanI 
ll992l : ISpringeill2010l ) code Gadget-2 (|Springelll2005l ). We 
first outline our method for setting up idealised disc gala:x- 
ies, before discussing the sub-grid physics used in the simu- 
lations. 



2.1 Initial conditions 

The full method for constructing equilibrium disc galaxies is 
given in AppendixiA here w e provide a summary. We largely 
fo UowlSpringel et al.l (|2005l ) which builds on the earlier work 
of lMo et al.l ( 19981 ). Our approach differs in a few key areas 
such as our treatment of adiabatic contraction of the dark 
matter (DM) halo (see Section IA1.3|I and the model used 
to describe the interstellar medium (ISM; Section [2]2]). Our 
model galaxies consist of a DM halo surrounding a disc made 
up of stars and gas, with a central stellar bulge. We assume 
that the DM halo and bulge follow Hernquist profiles and 
that the disc components are axisymmetric, following an 
exponential surface density profile with both gas and stellar 
components having the same scale length. The gas disc is 
assumed to be in hydrostatic equilibrium, whereas the stellar 
component is an isothermal sheet. 

Our galaxy parameters are as follows. The total mass 
is M200 = 1.36 X lO^^MfJB, with a Navarro, Frenk fc White 
(NFW: lNavarro et al.lll997l ) concentration parameter, c = 9. 
The mass is split (stated as fractions of M200) according to 
the disc fraction /d = 0.04 and bulge fraction /b = 0.01, 
with the remaining mass residing in the DM halo. (We do not 
include a gaseous halo in this study but note that one quickly 
forms when supernova feedback is included.) The disc is split 
once more into gas and stars with the mass fraction of the 
disc in gas, /g = 0.3. We set the radial disc scale length, 
h — 3.45 kpc, through angular momentum considerations 
and use it to define the bulge scale length and disc scale 
height b = zo — 0.2 h. Our galaxy parameters are chosen to 
approximately describe a Milky- Way -like galaxy and allow 
comparison with earlier work s uch as Springel et al. (|2005h : 
ISchave fc Dalla Vecchial l|2008l) and IDebuhr et al.l (|2010l '). 

The galaxies were constructed in both high and low 
resolution realisations and were run for 3Gyr. Table [T] sum- 
marises our choice of numerical parameters. The additional 
physical processes affecting the gas (over and above gravity 
and standard SPH) are discussed below. 



^ M200 is the mass enclosed within -R20O1 the radius at which the 
mean density is 200pcr, where per is the critical matter density 
required for a flat universe. 
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Table 1. Numerical parameters for our model galaxies, simulated at low and high resolution. Columns 2-6 give the total particle number, 
as well as the number in the DM halo, gaseous disc, stellar disc and stellar bulge respectively. Columns 7 & 8 list the dark matter and 
baryonic particle masses while Column 9 gives the (equivalent Plummer) gravitational softening length used in our simulations. All of 
our main results use the high resolution simulations. 



Resolution 


N 


A^h 


N^ 


JV* 


A^b 


moM [Mq] 


^baryon [Mq] 


Esoft [kpc] 


Low 
High 


7 X 10^ 
7 X 10"^ 


591089 
5910890 


21780 
217800 


60990 
609900 


26138 
261380 


3.08 X 10^ 
3.08 X 10^ 


6.24 X 10^ 
6.24 X 10"' 


0.05 
0.05 



2.2 Radiative cooling, ISM & star formation 

Gas is allowe d to cool radiatively following the method 
described by Kay 112004) . who adopted the approach of 
iThomas fc CouchmanI lll992|^. We utilise a tabulated equ ilib- 
rium cooling function from [Sutherland fc Dopital (|l993r ) for 
a metallicity Z — 0.3 Zq gas. The temperature of the gas is 
not permitted to fall below a critical value Tc — lO* K, at 
which point we assume the gas temperature is maintained 
by the meta-galactic UV background. 

The ISM and star formatio n model employed is that 
of ISchave fc Dalla Vecchial l|2008l '). We assume that the ISM 
is an ideal multi-phase gas which may be described by a 
polytropic equation of state above a critical value for the 
number density of hydrogen, nu,c 



Anr 



nnkBTc/iXfi) 



if riH > nH,c 
otherwise, 



(1) 



where ^ is a constant set by ensuring the pressure is con- 
tinuous across the boundary nn = wh.c and 7cfj — 4/3. The 
critical density is set to wh.c = O.lcm"^, the mean molecular 
weight, /i = 0.59 and the hydrogen mass fraction, X = 0.76, 
the latter two having primordial values. Gas may leave the 
equation of state when it receives (AGN or supernova) feed- 
back energy or if its t hermal energy increase s by 0.5 dex in a 
single timestep (as in lBooth fc Schavell2009 l. Allowing par- 
ticles to leave the equation of state in this fashion prevents 
artificial damping of feedback driven outflows. 

A gas particle lying on the equation of state may be 
converted to a coUisionless star particle representing a simple 
stellar population (SSP) described by a Salpeter initial mass 
function (IMF). The star formation ra te is designed to match 
the observed Kennicutt- Schmidt law (jKennicuttll 19981 ) for a 
galaxy with disc gas fraction, /g — 1. Specifically, the law 
used here is 



E* = 2.5 X lO"*M0yr"^kpc" 



IM0PC- 



(2) 



where E* is the star formation rate surface density and Eg is 
the gas surface density. Fig. [l] demonstrates the good match 
between simulated and observed star formation rates, once 
the latter is adjusted for the lower disc gas fraction (/g = 
0.3) adopted for our initial conditions Q 



^ We have also checked that our implementation reprodu c es the 
star formation rates found by ISchave &: Dalla Vecchial 1120081 ) 
when we adopt a zero metallicity gas, for both Salpeter and 
Chabrier IMFs. 




10" 10' 



Figure 1. Recovered Kennicutt-Schmidt relation for a high res- 
olution disc galaxy evolved for a time, t = 1 Gyr. Each point 
represents the star formation rate density averaged within an an- 
nulus centred on the galaxy rotation axis, while the colour repre- 
sents the age of these stars at the end of the simulation (with red 
points being the oldest stars). The dashed lines show the input 
Kennicutt-Schmidt relation and density cut-off, while dot-dashed 
lines show these adjusted for /g = 0.3 (the value used in our 
initial conditions). 



2.3 Supernova feedback 

The supernova feedback implementation described here is 
motivated by the desire to ensure high temperature gas 
is produced, as this is more likely to have an impact on 
the star formation rate. This is because gas with T ^ 
10^ K cools more slowly than cooler gas, mainly via ther- 
mal bremsstrahlung, allowing the heated material to in- 
teract hydrodynamically with the surrounding gas and es- 
cape the star forming region. We guarantee a high tempera- 
ture by reducing the mass loading factor of the heated gas, 
given a fi xed amount of energy availab le from the super- 
novae (see iDalla Vecchia fc Schave 2012 for an application 
of this method to disc galaxies and lKav et al. 2003 for galaxy 
groups) . 

We assume that a fraction of the stars in each SSP 
will end their lives as type II supernovae (SNII), releasing 
energy to the environment. This typically occurs on a short 
timescale which we assume to be negligible. The amount of 
energy provided by SNII per unit mass of stars formed is 
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calculated as 

USN = EsN nsNII, 

= -EsN / $(M) dM, 

J Mo 



from a surrounding gas reservoir, given by 



5.6 X lO^^ergg"^ 



-EsN 
IQSierg 



nsNii 



1.11 X 10-2Mq^ 



(3) 



where Esn is the energy per Type II supernova and nsNii 
is the number of supernovae per unit mass, obtained by 
integrating over the IMF ($) between Mo — 0.1 Mq and 
Ml = 100 Mq . The specific energy change imposed on the 
local gas is then 



Ithoat ~ eSN 



rrista 



?TlgasA*'hcat 



USN, 



(4) 



where esN is a dimensionless supernova efficiency parameter 
(accounting for small-scale energy losses) and the term in 
brackets is the inverse mass loading factor (heating A'hcat 
gas particles of mass TTigas for every new star particle with 
mass TTistar). The heating temperature can then be expressed 
as 



Theat ~2.7x lO'K 

-Esn 



esN 

A^hoat 



V"^ga: 



V 1051 erg 



0.59 



nsNii 



1.11 X 10-2Mq^ 



(5) 



For the purpose of our fiducial simulations we set Thcat = 
lO'^ K and TVhoat = 1, yielding an assumed efficiency of 
egN ~ 0.4 (this choice is investigated further in Section [5. 3p . 
When heated by supernovae, a gas particle is removed from 
the equation of state on its current timestep and tagged as 
a wind particle. However, should the particle meet the con- 
ditions for joining the equation of state on any subsequent 
timestep, it will return as normal. 



2.4 Black hole growth and AGN feedback 

There exists a large variety of work employing different 
methods for simulating AGN and their effects on the host 
galaxy. In this paper we study three of the most wid e ly use d 
imple mentations, namely those by ISpringel et al.l l|2005f ): 
iBooth fc Schaye (2009i ) ; and a novel method based on the 
model of iPower et alT l|201ll 'l. The AGN models may be 
simplistically broken down into three processes: formation, 
growth and feedback. The formation of black holes at high 
redshift is an interesting topic in itself and an active area 
of research. However, it is beyond the scope of this pa- 
per where we simply assume our galaxies already contain 
a super-massive BH at their centre. We discuss the other 
two processes in turn. 



2.^.1 Black hole growth: the Bondi method 

ISpringel et all (|2005l ) and iBooth fc Schavd (|2009l ) both 
used modified versions of the B ondi-Hoyle-Lyttleton 
JHovle fc Lvttletonll 19391 : iBondi fc Hovle 1944 : Bondi 1952 ) 
model to calculate the accretion rate onto the black hole 



Mbh = a 



47rG^MgHpgae 

(c2+„2)3/2 



(6) 



where pgaa is the density of the surrounding gas, Cs its sound 
speed, and v the velocity of the B H relative to th e gas. T he 
additional factor a was added by ISpringel et al.l (12005! ) to 
counteract an assumed underestimate for the accretion rate 
due to inadequate numerical resolution. Simulating the ac- 
cretion fully requires proper resolution of the Bondi radius 



GMb 



re 



(7) 



which is considerably smaller than the force resolution in 
current simulations, for the case of galactic-scale black holes. 
Furthermore, one must also resolve the multiphase ISM and 
crucially its cold dense clouds, to allow accurate computa- 
tion of_th£_gas_d^nsity and sound speed close to the BH. 

[Booth fc Schavd (l2009J) modified the correction factor, 
a, arguing that the accretion is not always equally under- 
resolved and suggested the following scaling with density 



if TlH < nH,c 

otherwise. 



(8) 



where tih.c is the same critical value as for the ISM and 
star formation (jih.c = 0.1 cm~^). The i mpact of varying 
the po wer-law index, /?, was investigated bv lBooth fc Schavd 
(J2009f ) as well as varying the constant a value, as used in 
the more standard formalism. They found that observables 
are much more sensitive to the choice of a than to any rea- 
sonable change in /3. We have implemented both methods 
as laid out in their original form. 

A BH particle's mass is tracked by two variables: the 
dynamical mass which increases instantaneously when a par- 
ticle is captured by the BH; and the smoothly integrated in- 
ternal mass which is used to determine the accretion rate of 
the BH. In practice, the BH particle grows by stochastically 
capturing surrounding gas particles within the smoothing 
length at a rate which coarsely ma tches the internal black 
hole mass (see e.g.lDi Matteo et al . 2005; Sii acki et "aLll2007l : 
iDi Matteo et al.ll2008l : iBooth fc Schavc.,20091 ). Upon accret- 
ing a particle the mass is added to that of the black hole. The 
potential drawback of using the SPH kernel is that gas parti- 
cles may be captured at arbitrarily large distances from the 
sink particle due t o the use of an adaptive smoothing len gth 
(discussed e.g. bv lDebuhr et al.ll20ld : IPower et aLlbOllI ). 



2.4-2 Black hole growth: the accretion disc method 

The IPower et al.l (|201lj) method does not use the Bondi- 
Hoyle-Lyttleton equation. Instead, it assumes that any par- 
ticle which passes within a fixed accretion radius, race, is 
captured onto the accretion disc surrounding the BH. This 
mass is then moved to the BH on a viscous timescale. A 
key feature of this method is that, unlike the Bondi-Hoyle- 
Lyttleton accretion, only gas with low angular momentum is 
likely to be captured. In simulations lacking sufficient spatial 
resolution to s imulate true accr etion discs existing on sub- 
pc scales (e.g. iKato et al.lll998l ). it is favourable to set race 
to the smallest resolvable length-scale, i.e. the gravitational 
softening length Esoft. 
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Table 2. Summary of our main runs and important parameter choices. Column 1 lists the label used to refer to each model, while 
columns 2 & 3 specify whether supernovae and AGN are included. Pertinent details of the run are given in column 4; columns 5-7 lists 
values of the appropriate BH accretion rate parameter; column 8 the number of particles heated by the BH in one event and column 9 
the assumed value for the feedback coupling efficiency parameter. The value of A^hcat fo'' the KAGN runs is approximate. 



Label 



SN AGN Details 



/3 



fff 



N^ 



NFB 


N 


N 


SN 


Y 


N 


KAGN 


N 


Y 


KAGN.SN 


Y 


Y 


SAGN 


N 


Y 


SAGN.SN 


Y 


Y 


DAGN 


N 


Y 


DAGN.SN 


Y 


Y 



No feedback 

Supernovae only 

Bondi accretion, kernel feedback 

Bondi accretion, kernel feedback, supernovae 

Boosted Bondi accretion, strong feedback 

Boosted Bondi accretion, strong feedback, supernovae 

Disc accretion, strong feedback 

Disc accretion, strong feedback, supernovae 



100 


- 


- 


50 


0.05 


100 


- 


- 


50 


0.05 


- 


2 


- 


10 


0.15 


- 


2 


- 


10 


0.15 


- 


- 


0.1 


10 


0.15 


- 


- 


0.1 


10 


0.15 



Preliminary testing of this method found it to be un- 
physical at the typical resolution of our simulations (~ 
100 pc), yielding an artificially large accretion rate, an 
overly-massive black hole and a large clear e d-out a rea at the 
centre of the galaxy. IWurster fc Thackeij l|2013bl ) attempt 
to overcome this problem by significantly reducing race to 
limit mass accretion, with the precise value being varied as 
part of their investigation. It is, however, unphysical to take 
particles as point masses in this regime since capturing one 
whole corresponds to a gas cloud with radius of order the 
SPH smoothing length instantaneously collapsing onto the 
much smaller accretion disc. We take an alternative route 
by instead interpreting the accretion radius as a scale within 
which the captured material is gravitationally bound to the 
centre of the galaxy and will contribute to the mass of the 
accretion disc on a timescale set by the local free-fall time, 
iff = l/\/G(pgas). Specifically, the gas mass is added to the 
disc at a rate given by 

dMdisc 



effMgas(< race) 



(9) 



dt iff 

where eg is an efficiency parameter that may represent un- 
resolved effects such as turbulence or residual angular mo- 
mentum (we set egf = 0.1, as discussed in the next section). 
Motivating the larg e scale accretion via the freefall timescale 
is similar in spirit to lHobbs et al.l (|2012h , who did so by mod- 
ifying the Bondi model. Mass is then added to the BH at 
the rate 



Mbh = rnin 



Mdi„ 



i-vis' 



,AfE 



(10) 



where the viscous timescale fvisc rnay be estimated from 
physical arguments and is typically set to fvisc ~ 10 — 
100 Myr (JKato et al.ll 19981 : IPower et alll201ll : we set fvisc = 
10 Myr). For all three methods, we restrict the accretion rate 
to the Eddington limit 

• inGMBHmH 

JWEdd = , (11) 

CrCTTC 

where we assume the standard value for the radiative effi- 
ciency parameter, er =0.1. 

Independent of accretion r ate method, particles are nu- 
merically captured similarly to lSpringel et al.l (|2005h . Parti- 
cles are accreted only when the system mass (BH and accre- 
tion disc if applicable) exceeds the BH particle's dynamical 
mass. Accreted particles would ideally contribute momen- 
tum to the BH particle as it is clearly desirable to conserve 



momentum in all simulations. However we found that arti- 
ficially large kicks and two-body scattering from other par- 
ticles caused the black hole to become displaced from the 
centre of the galaxy, even in iso lated quiescent galaxies . We 
therefore apply the app roach of lDi Matteo et al.l (|2005J) and 
iBooth fc Schavd l|2009r ) to all three methods, and re-centre 
the black hole particle on the neighbouring (r^j < h) particle 
of minimum potential until it is a factor of 10 more massive 
than a single gas particle. 



2.4.3 AGN feedback 

The accretion rate onto the central BH is used to determine 
the energy budget for the AGN feedback. It is commonly 
assumed that a fraction of the accreted rest mass energy is 
released from the BH, with an additional multiplicative fac- 
tor to represent the efficiency with which this then couples 
to the surrounding gas 



Ef, 



— efLr = efErMBHC 



(12) 



where e^ = 0.1 is the radiative efficiency as before and ef is 
the feedback coupling efficiency. 

The most accurate way to distribute this feedback en- 
ergy, and the form it takes, is unclear. One method is to 
utilise an SPH-like approach and simply deposit the energy 
thermally in a kernel-weighted fa shion; such an approach 
was taken bv lSpringel et al.l (|2005J). This method is numeri- 
cally simple and characterises a high energy quasar radiation 
field coupling isotropically to the surrounding gas. However, 
releasing the feedback energy in this way may not be able to 
provide gas hot enough to dri ve strong outflows . An a lter- 
native approach was taken bv iBooth fc Schavd (|2009r ). By 
'bottling up' the energy it is ensured that a critical tem- 
perature change, ATmin, is reached before a heating event 
occurs. The energy which must be stored to heat A'hcat gas 
particles by ATmin is 

„ A'"hcat»ngas fee ATmin ,.,„-, 

iicrit = -. -T • (13) 

(7 — l)/imH 

Once heated, it is then assumed that the gas will rise buoy- 
antly in the ambient medium, mimicking outflows found in 
observations. 

In the following sections we perform simulations that 
employ the methods for describing gas processes and feed- 
back outlined above. Specifically, we simulate the Bondi ac- 
cretion method with kernel-weighted feedback as set out by 
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log(S/M„„kpc-^) 

Figure 2. Projected face- and edge-on gas surface density maps for an isolated galaxy simulated with various feedback implementations, 
at t = 0.6 and 3Gyr. Columns show results for the NFB, SN and SAGN models respectively. The panels are 100 kpc on a side and 
surface density is coloured on a log scale. The reduction in density with time, seen for all models, occurs due to star formation. A slight 
redistribution of gas is seen for the SAGN model but the most pronounced case is for supernova feedback (SN), where a hot halo is 
quickly generated. 



ISpringel et"ai1 (|2005l ). labelled he re as KAGN; the Bond i 
method with 'strong' feedback by iBooth &: Schavd l|2009l ). 
labelled SAGNj_and the disc accretion method based on 



IPower et al] (|201ll ') with 'strong' feedback (DAGN). We also 
consider a model with no feedback (NFB) and a model with 
supernovae only (SN). Models run with both AGN and su- 
pernovae are denoted XAGN_SN, where X is (K,S,D) as dis- 
cussed. Table[2]summarises the details of these runs and lists 
values for the main AGN feedback parameters. 



tions of an isolated disc galaxy provides us the opportunity 
to study the models in a regime where the black hole is 
typically expected to be in a phase of low activity and pro- 
vides a comparison for the behaviour of black holes in violent 
mergers. Employing such high resolution 'toy models' also 
ensures any observed phenomena may be more easily dis- 
entangled because the systems themselves are simple and 
relatively well understood. 



3 ISOLATED DISC GALAXIES 

We begin our investigation by performing tests of the feed- 
back recipes in isolated disc galaxies. Performing simula- 



3.1 Simulations with and without supernova 
feedback 

We first investigate the effect of supernovae on the host 
galaxy properties, neglecting the AGN for the time being. 
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Fig. [2] shows projected surface density maps for the NFB 
and SN runs at i = 0.6 and 3Gyr. The maps illustrate 
qualitatively that supernovae quickly generate a gas halo 
that persists for the duration of the simulation. They also 
modify the disc structure, increasing the gas porosity. The 
reduction in density in the disc over time is predominantly 
due to star formation as this can be seen for the NFB case. 

Fig. |3] shows the evolution of the global star formation 
rate for the NFB and SN simulations (solid and dashed lines 
respectively). We also show the star formation suppression 
factor, Sx = M*(NFB)/Af*(SN), in the bottom panel. As 
expected, we see a reduction in the star formation rate when 
supernovae are included. This persists until gas depletion 
becomes important (i ~ 1.5 Gyr) when we see an increase 
in star formation. This is due to low-level star formation 
from gas which would have otherwise already formed stars 
in the absence of supernovaqj. 

The primary change seen in simulations with supernova 
feedback is a reduction in the stellar mass of the disc. This is 
clearly seen in Fig. |31 where both the mass in new stars and 
in a supernova-driven wind is shown as a fraction of initial 
gas mass, versus time for the NFB and SN simulations. A 
reduction of ~ 10 per cent in stellar mass is seen by the end 
of the simulation, caused by the inclusion of supernovae. It is 
clear that the majority of supernova heated gas is returning 
to the equation of state on a short timescale, as the mass in 
the wind is considerably lower than in new stars. 



3.2 Simulations with AGN feedback only 

We now consider runs with AGN feedback only. Unlike the 
supernovae, the AGN feedback has no effect on the star for- 
mation rate for any model (an example is shown in the bot- 
tom panel of Fig. [3] for the SAGN run). This is likely due to 
the differing scales involved as the black hole typically only 
interacts with the central ~ 0.1 kpc whereas the star for- 
mation takes places on kpc scales. Although the black holes 
do not change the large-scale galaxy properties they can be 
seen to drive a weak wind. 

The evolution of BH mass with time is shown in the top 
panel of Fig. O where we can see that the three models each 
produce very different final masses. However, we note that 
the normalisation and therefore the final mass may easily 
be tuned through accretion rate parameter choices for any 
given model. We have chosen to employ the original values 
where appropriate (KAGN and SAGN) as these were tuned 
to recover the Mbh — Afbuige re l ation in full cosmological 
simulations. For the lPower et al.l (|201ll ) method (which was 
not tested on cosmological volumes) we adopted the orig- 
inal value for ivisc(= lOMyr) and tuned the free-fall efS- 
ciency parameter, eg = 0.1, so as to bri ng the final mass 
roughly in line with the lBooth fc Schavd (2009) method as 
we adopted their feedback prescription for the DAGN model. 
The most significant factor affecting the final BH mass in 
isolated galaxy simulations is the resolution used (we inves- 
tigate this in Section [53}. 



^ The suppression in t he star formation rate is found to be con- 
sistent with that seen in lDalla Vecchia fc Schavj 1(2013) under the 
appropriate parameter choices (see their Fig. 10). 
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Figure 3. Star formation rate (top) and suppression factor (bot- 
tom) throughout simulations of an isolated high-resolution disc 
galaxy. Suppression is calculated as a ratio of instantaneous star 
formation rates, no-feedback case over feedback case. Supernova 
feedback causes the dominant suppression of the star formation 
rate, while the AGN wind is unable to interact with the majority 
of the galaxy, only having a minor influence when able to couple 
with the supernova-inflated hot halo. 



We also note that the models predict different BH his- 
tories. Bondi-based methods predict that the BH grows 
smoothly whereas in the DAGN method it undergoes an 
initial spike with an abrupt transition to a period of slower 
growth. Unlike the spatially adaptive Bondi methods, the 
disc method employs a small and finite accretion radius 
which initially contains gas (an artefact of our initial con- 
ditions) but is quickly removed through accretion and feed- 
back, shutting off the black hole growth. After this time, the 
BH growth rates are similar in the SAGN and DAGN mod- 
els, while the BH in the KAGN method continues to grow 
at a faster rate. This suggests that the choice of feedback 
plays an important role in determining the BH growth rate 
(we confirm this in Section [5. 6|l . 

It is also interesting to note that none of the black holes 
in these simulations reach the mass scales (^ 10^ — 10^ IVIq 
predicted from the ob servations (e.g. Magorrian et al.lll99i 
iTremaine et al.ll2002l ). There are some possible causes for 
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Figure 4. Fractional mass in new stars and particles tagged 
as being in a supernova wind versus time for the isolated high- 
resolution disc galaxy, with and without supernova feedback. All 
masses are shown as a fraction of the initial gas mass. Super- 
novae act to reduce the amount of gas on the equation of state 
and therefore the mass in stars. 
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Figure 5. Evolution of BH mass with time, for the fiducial AGN 
models without supernovae in isolated galaxy simulations (top 
panel). The ratios of masses found between these and correspond- 
ing simulations including supernovae are shown in the bottom 
panel. Supernovae increase the BH mass in all models by pre- 
venting gas depletion through star formation and driving gas to- 
wards the BH. However, this increase (~ 10 per cent) is much 
smaller than the differences between the AGN feedback models 
themselves (which can be tuned if necessary). 



this discrepancy, several of which are discussed later in this 
paper, such as reduced accretion in the absence of mergers, 
low initial black hole mass, resolution effects and large-scale 
gas accretion. Investigating th e final issue would require the 
inclusion of a hot gas halo (e.g. Sinha fc HoUev-Bockelmannl 



[2009"; 'Mos ter et a"l1l201ll : IWurster fc Thackeijl2013bl ) and/or 
the simulation of l arge-scale accretion onto the galaxy (e.g 
iMoster et al.ll2012l ). and is beyond the scope of this paper. 



3.3 Simulations with supernova & AGN feedback 

We now investigate the effect of combining both supernova 
and AGN feedback in our simulations. Fig. |3] shows a mild 
reduction in the star formation rate in the simulations in- 
cluding both strong AGN and supernovae (SAGN_SN), com- 
pared with the supernova only case (SN). This is due to AGN 
winds coupling to a supernova-fuelled hot gas halo reducing 
cooling and subsequent star formation (the effect is slightly 
more prominent in lower resolution simulations; not shown). 
This has been verified by examining the least dense parti- 
cles in the simulations and it was found that gas driven out 
to large radii by strong AGN feedback is missing in simu- 
lations also including supernova feedback. The effect is not 
seen in the KAGN_SN models as there the AGN drives a 
much weaker wind due to the larger mass being heated. 

Semi-analytic galaxy formation models commonly as- 
sume that the effects of supernova and AGN feedback af- 
fect galaxy properties independently of one anothe r (e.g. 
iBower et al.ll2006l : ICroton et al.ll20od : lGuo et al.|[201ll ). This 
approach has recently been brought into question by the 
work of Booth & Schayc (2012), who found that the two 
feedback processes suppressed the star formation rate less 
than if they were assumed independent. For a more quan- 
titative analysis we ca lculated the feedback su ppression 'x' 
parameter as used in iBooth fc Schavel (|2012l ). where x is 
defined as 



X- 



Sag 



N+SN 



Sagn X 5's^ 



(14) 



and Sx is the suppression in the star formation rate due to 
process X. The value of x quantifies the level of interaction 
of AGN and supernova feedback which results in the sup- 
pression of star formation. A value of x > 1 (x < 1) shows 
that the two feedback processes have amplified (weakened) 
each others' ability to reduce the star formation rate and 
X = 1 indicates that the two processes are independent. We 
find X ~ 1 at all times for KAGN runs and x — 1 rising 
to X - 1-15 at late times for the SAGN and DAGN runs. 
This indicates that the interaction between the two feedback 
processes is of minimal importance, even for the strong feed- 
back simulations where a coupling is observed. This weak 
amplification disagrees with the value of x ~ 0.3 found by 
iBooth fc Schavel (|2012l ) for systems of M ~ 10^^ Mq such 
as this. However as this work considers one high resolution 
object in detail whereas iBooth fc Schavel (|2012l ) perform a 
cosmological simulation, this finding is not necessarily a di- 
rect contradiction. Both the environment of the gala^cy and 
how well the outflow is resolved are likely to be important 
in determining the strength of this effect. 

The lower panel of Fig. [S] shows the ratio of BH masses 
for the main AGN models with and without supernovae over 
the course of the simulations. Regardless of model, the black 
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Figure 6. Projected gas surface density maps for merging galaxies in the NFB (top panels) and DAGN simulations (bottom panels). 
Columns show snapshots taken at times, i = 0.3,0.6 and 3 Gyr, viewed perpendicular to the merger plane. Panels are 250 kpc on a side 
and colours represent surface density on a log scale in units Mq kpc~^. Viewed at i = 0.3 Gyr, the upper (lower) galaxy has entered 
from the right (left) with its rotation axis tilted (parallel) with respect to the line-of-sight and swings around to appear as the lower 
(upper) galaxy at t = 0.6 Gyr producing tidal tails and bridges, before finally coalescing. 



hole mass is mildly enhanced by the action of supernovae 
after 1 Gyr due to a combination of reduced star formation 
leaving more gas to accrete and winds from the star-forming 
disc feeding the central BH. 

In summary, the reduction in star formation rate in our 
isolated galaxy is dominated by the effect of supernovae, 
whilst AGN have very little impact regardless of whether the 
feedback processes are simulated in isolation or in tandem. 
The final black hole mass is found to differ between models 
and consistently shows a mild increase as a result of the 
inclusion of supernovae. 



4 MERGING DISC GALAXIES 

Having investigated the influence of supernova and AGN 
feedback processes on isolated systems, we now turn our at- 
tention to how they interact through a major merger. It is 
important to consider merging systems as they are thought 
to be a critical phase in galaxy evolution, driving both peri- 
ods of high A GN activity and starbursts. It has been shown 
previously by ISpringel et al.l (|2005|) that the inclusion of 
AGN feedback terminates star formation and halts the star- 
burst. Here we re-investigate this result to establish whether 
this is a generic feature seen in our runs. 



4.1 The merger scenario 

We have chosen to simulate a 1:1 major merger between 
identical disc galaxies as described in Section [2. II This ex- 
treme choice was made primarily to allow for easier com- 
parison with the isolated galaxy simulations and with pre- 
vious work. We start the merger with an initial separation 
of 150 kpc, and a perpendicular impact parameter of 25 kpc. 
The galaxies initially have approximately zero total energy, 
corresponding to a net velocity for each galaxy of 185 kms~^ 
in the centre of mass frame, directed parallel to the separa- 
tion axis. One of the galaxies is rotated by 30° with respect 
to both the axis of symmetry and the plane of the merger. 
This merger geometry is i ntentionally similar to earlier work 
fe.g. lSpringel et al.ll2005l : iDebuhr et al.ll2010l ). Our choice of 
merger set-up results in the galaxy centres experiencing two 
passes, coalescing on the second, and merging at around 
t ~ 1.5 Gyr. Three key phases of the merger are illustrated 
in Fig. [G] which shows the column density of gas for the first 
pass, largest subsequent separation and final state, for the 
NFB and DAGN simulations. 

At i = 3 Gyr the merging system has had time to relax 
back into an equilibrium state. The stellar population of the 
remnant is a slowly rotating oblate ellipsoid with sphericity, 
s(= c/a) = 0.49; elongation, e(= b/a) = 0.86; and triaxial- 
ity, T[= (a^ - 6^)/(a^ - c^)] = 0.33 where o > 6 > c are the 
square roots of diagonal elements in the mass distribution 
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Figure 7. Star formation rates and suppression factors for merging galaxy simulations with supernova feedback (left panels); AGN 
feedback (middle panels); and both (right panels). Vertical lines indicate the times of first and second passage. The suppression factor is 
the ratio of the feedback case to the no-feedback case. Strong AGN feedback has the largest impact on the star formation rate, eliminating 
the main star burst. 



tensor (see iBrvan et al.ll2012l ) . The abundance of gas at late 
times is model-dependent, but takes the form of a slowly 
rotating extended halo with a rotating {vc ~ 200kms~^) 
dense gas disc at its very centre. The central disc is formed 
at a small angle to the plane of the merger due to the initial 
tilt applied to one of the galaxies (Fig. |6]|. All merger rem- 
nants have low stellar angular momentum and are classified 
as 'slow rotators' under the sche me appli ed to the ATL AS 
3D galaxies (e.g. lEmsellem et al.i 2007 : Bois et al.ll201ll) in- 
dependent of any feedback processes. 

4.2 Mergers with supernova feedback only 

We begin, as with our investigation into isolated galaxies, 
by evaluating the effect of the supernovae. The evolution of 
the star formation rates for the main simulations are shown 
in Fig. [7] as well as the suppression factor relative to the 
no-feedback case. For the NFB and SN models, an increase 
in the rate is seen for both passes (t ~ 0.3 and 1.3 Gyr), 
however the burst at the time of second passage is substan- 
tially more significant, as found in previous studies when 
the galaxies contain a bulge (e.g. iMihos fc HernauistlllQQB : 
ISprineel et al.ll2005l : lDebuhr et al.ll2010h . The primary burst 
takes place just as the galaxies coalesce, when the gas discs 
collide most strongly. As in the isolated case, we find that 
supernovae reduce the star formation rate throughout the 
majority of the simulation. However, star formation through 
the starburst is unchanged, suggesting that the supernova- 
heated gas is largely confined during this period. 

Fig.[8]shows the fraction of initial gas mass turned into 
stars versus time for the NFB and SN simulations. We also 
plot the fraction heated by supernovae and gas currently on 
the equation of state for the latter case. Almost all of the 
gas is converted into stars by the end of the simulation in 
the NFB run while supernovae reduce this fraction by ~ 13 
per cent, similar to that seen in the isolated galaxy simula- 
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Figure 8. Fractional mass in new stars and gas in various phases 
for simulations of high-resolution disc galaxy mergers with and 
without supernova feedback. For the NFB case only the stellar 
mass is shown, while for SN we additionally plot the supernova- 
heated particles and gas on the equation of state. All masses are 
shown as a fraction of the initial gas mass. As in the isolated case, 
supernovae reduce the amount of gas on the equation of state and 
stellar mass formed. 



tions. Supernovae cause an additional reduction of mass on 
the equation of state early on, reducing star formation and 
preserving gas. The amount of gas on the equation of state 
temporarily increases twice, coinciding with the timing of 
the galaxy passages that cause an increase in gas density 
during the collision. This leads to the enhanced star for- 
mation rate that acts to reduce the amount of gas on the 
equation of state once more. The amount of gas previously 
heated by supernovae shows only a weak change at the times 
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Figure 9. Final temperature distributions versus radius for gas particles in the NFB, KAGN_SN and SAGN_SN (left to right) simulations, 
calculated at t = 3 Gyr. Symbols indicate the 'tags' applied to gas particles. The predominant trend in all panels is the same, indicating 
that the dynamical effects of the merger cause much of the heating and redistribution of the gas. However, the runs with strong AGN 
feedback (SAGN_SN is shown here) produce a larger amount of hot gas (T > 10^ K) at large radii (r > 100 kpc). 



of the merger passages as it has escaped to areas of lower 
density. The number of particles tagged as heated is, as in 
the isolated simulation, only a small fraction of the num- 
ber of stars formed because the majority quickly rejoins the 
equation of state without escaping the disc. 

We present in Fig. [9l the final temperature distribu- 
tions of gas particles for the NFB, KAGN^N and SAGN_SN 
simulations (left to right). Here we discuss the effect of su- 
pernovae using the KAGN_SN particle temperature distri- 
bution as a proxy for SN, to avoid duplication as they are 
remarkably alike. The distributions for NFB and SN are 
very similar, producing a halo of diffuse, hot gas that ex- 
tends to ~ 1 Mpc from the halo centre. This indicates that 
the dominant cause of heating is gravitational in origin and 
a consequence of the merger. The importance of dynamical 
effects in the formation of the hot halo is corroborated by 
the increase in gas fraction in the halo. This is 53 per cent 
at the end of the simulation when feedback is omitted com- 
pared to 0.02 per cent for the isolated galaxy (fractions with 
supernovae increase to 90 and 48 per cent respectively). 

The final gas density profiles are shown in Fig. 1101 Pro- 
files are shown at t = 3 Gyr and centred on the centre of 
stellar mass; also shown are the ratios of the density pro- 
files compared to the NFB simulation. The SN profile shows 
only a small increase in gas density beyond 10 kpc due to a 
reduction in star formation. 



4.3 Mergers ^vith AGN feedback only 

In stark contrast to the isolated galaxy case we find that in 
merging systems, the AGN can have the dominant effect on 
the star formation rate when the feedback is strong (SAGN 
and DAGN), even eliminating the main starburst (Fig. [T]). 
Both methods greatly reduce star formation whereas kernel- 
weighted feedback (KAGN) does not. However the star for- 
mation rate appears to be insensitive to the choice of BH 
accretion rate (again comparing SAGN and DAGN). Over 
the whole simulation, the weaker kernel feedback (KAGN) 
reduces the final mass in stars by only 5 per cent, by reduc- 
ing the starburst duration and strong feedback by ~ 25 per 
cent. We investigate this further in Section [5.61 

The wide disparity in the power of the AGN to infiu- 



ence its environment is likely due to the differing temper- 
atures which gas is heated to. At 10* K (the strong feed- 
back heating temperature) radiative cooling has less of an 
effect, as the cooling time increases as -v/T when the emis- 
sion is dominated by thermal Brehmsstrahlung. Therefore, 
by spreading the available feedback energy over many parti- 
cles neighbouring the black hole and depositing any available 
energy promptly, the kernel-weighted method heats more 
gas to lower temperatures where it can cool more easily.The 
centre and right panels of Fig. [9] show that whilst kernel 
feedback removes gas from the equation of state and up to 
temperatures of ~ 10^ K, strong feedback heats much more 
gas to higher temperatures and larger radii (up to 10® K at 
10'' kpc), resulting in less gas on the equation of state. 

Final gas density profiles (middle panels of Fig. llOfl 
show a redistribution of gas out to large radii in the case of 
the strong feedback models. The kernel-weighted feedback 
method makes very little difference to the gas distribution, 
causing only a slight increase in gas at r < 10 kpc. The 
strong feedback methods however redistribute gas out to 
much larger radii, r > 10^ kpc, with an order of magnitude 
larger density than stars at that radius, increasing the mass 
of gas at r > 10'^ kpc by two orders of magnitude from 2 x 10^ 
to 2.5 X 10^ M0. This does not, however, significantly change 
the baryon fraction within -R200 compared to the NFB case. 
AGN have little effect on stellar density profiles, as found for 
isolated galaxies, due to the subdominant mass in new stars 
and the large amount of stellar mixing which occurs though 
a merger. Furthermore, gas removed to lower densities by 
strong AGN feedback is less likely to return and will never 
form stars over the course of the simulation. 

Evolution of the black hole masses with time for the 
three main AGN models are shown in Fig. 1111 Each panel 
shows the masses of the two black holes in the simulation 
as well as their sum. The black holes undergo a spike of 
growth at the times of the first and second passes, forming 
a tight binary system during the latter which has not coa- 
lesced by the end of the simulation (black holes are found 
to merge more readily in low resolution simulations). Cor- 
respondingly there is a large increase in the net amount of 
feedback energy deposited compared to the isolated case, 
which in itself suggests the AGN will influence its environ- 
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Figure 10. Final gas density profiles for the main set of simulations, grouped into runs with no AGN (left panels), AGN only (centre 
panels) and supernovae plus AGN (right panels). The lower sub-panels show the ratio for each model relative to the NFB run in the left 
and centre panels, and the SN run in the right panel. Both supernovae and strong AGN feedback affect the gas distribution, with the 
latter pushing gas out to larger radii. 



ment more strongly. The AGN activity peaks are coinci- 
dental -with merger passes, as the latter drives an inflo-w of 
gas, increasing the density around the black hole, and thus 
boosting the accretion rate. 



4.4 Mergers with supernova &: AGN feedback 

The previous sections have shown that AGN can play a more 
significant role than supernovae in determining the star for- 
mation rate during a merger. We can see from the star for- 
mation rates and ratios in Fig. [7] that AGN feedback -with 
supernovae is even more effective at reducing star forma- 
tion for the majority of the simulation. This is due to high 
temperature AGN-heated gas coupling to the gas halo in- 
flated by supernovae, slowing cooling and subsequent star 
formation. 

For a more quantifiable measure of the combined sup- 
pression effects of the feedback processes, we again calculate 
the X parameter. We find, as in the isolated galaxy simula- 
tions, X values consistent with unity prior to coalescence 
but rising to ~ 3 (1.4) at late times for the strong (kernel- 
weighted) models respectively. This, as in t he case of isolated 
galaxi es, is a departure from the findings of lBooth fc Schavd 
l|2012l ). However, we note that x values calculated after the 
merger are prone to noise due to the low-level of star forma- 
tion at this time. 

The coupling of the AGN outflow to the supernova- 
fueled hot gas halo, observed in both our isolated and 
merging galaxy simulations, raises the question of whether 
such a halo ought to have been included in our initial 
conditions. Such an approach has been adopted by recent 
galaxy merger studies (e.g. ISinha fc Hollev-BockelmannI 
12003 : iMoster et al.ll2011^ and may have affected our findings 
through the infiuence of feedback on late-time gas accretion 
from the halo. However, since we compare our galaxy models 
primarily to each other rather than observations any effect 



should be reduced, and simulations including supernovae 
quickly generate a halo in any case. Furthermore, omitting 
this component has allowed a high resolution study on the 
role of feedback processes which is more easily compared 
with preceding studies. 

Supernovae have a stronger effect on black hole mass in 
mergers compared to the isolated galaxy simulations. In the 
KAGN^N and DAGN_SN simulations the reduction in star 
formation resulting from the inclusion of supernovae results 
in more gas in the galaxy after coalescence to feed the black 
hole (Figs. [TOl and [TT|) . The heating from supernovae there- 
fore acts to increase growth by reducing star formation and 
leaving more surrounding gas at late times to be accreted 
by the black hole. This trend is not seen in the SAGN_SN 
simulation because the black holes merge promptly at coa- 
lescence, resulting in a comparatively reduced total accretion 
rate as the merged black hole stablilises at the centre of the 
system. 

In this section we have shown that, in contrast to simu- 
lations of isolated galaxies, it is possible for AGN to play the 
dominant role in suppressing star formation in a merger sce- 
nario. The method by which feedback energy is deposited is 
crucial in determining the strength of suppression; whilst 
models incorporating strong feedback eliminate the star- 
burst, the weaker kernel- weighted mode merely reduces its 
duration. Supernovae still play a role in suppressing the star 
formation during the majority of the simulation, but are 
found to be incapable of reducing the strength of the star- 
burst by themselves. 



5 A CLOSER LOOK AT THE MODELS 

Our main simulations have allowed us to analyse the role 
feedback processes play in disc galaxies for a range of exist- 
ing models. We now go on to look into numerical effects and 



© 2013 RAS, MNRAS OOO.IUMI 



AGN and supernova feedback in simulations of disc galaxies 13 




KAGN SN/KAGN 



SAGiN_SN/SAGN 



DAQN_SN/DAGN 




1.5 

f |Gyr] 



1.5 

t [Gyrl 



0.5 1.0 1.5 2.0 2.5 

t[Gyrl 



Figure 11. Evolution of black hole mass through a merger for the fiducial AGN models; left to right, KAGN, SAGN and DAGN 
respectively. Top panels show simulations without supernovae, where dotted and dashed lines show the black hole 'internal' smoothly 
integrated masses for the two black holes in the simulation and the solid line shows the sum of the two. Lower panels show the total 
black hole mass for each AGN model with supernovae relative to the no-supernovae equivalent. The horizontal black dotted line in the 
top panel shows the threshold mass above which black hole recentring is halted. Boosted growth rates can be seen at the times of first 
passage and coalescence of the merger. The inclusion of supernovae increases the final black hole mass when the two black holes do not 
merge (as is the case with KAGN and DAGN). 



the impact of our chosen parameters. As the high resolution 
achieved in these simulations is currently out of reach of 
typical cosmological simulations, we begin by investigating 
the impact of employing a lower mass resolution in Section 
15.11 The effect of increasing the galaxy gas fraction to bet- 
ter represent the high redshift systems in which most quasar 
growth and mergers are believed to occur is examined in Sec- 
tion 15.21 To analyse the consequence of choosing an initial 
black hole mass of 10^ Mq , we perform tests in which we set 
this to a larger value in Section [5.31 

We also discuss issues raised whilst performing this 
work and further investigate aspects of the feedback mod- 
els. Simulations of isolated galaxies showed that black holes 
appear to self regulat e, slowing their growth as found by 
IDj Matteo et al.l ( 20051 ) . However this may equally be caused 
by star formation consuming the gas which would otherwise 
feed the black holes; we study this further by artificially 
turning off star formation in Section 15.41 A further investi- 
gation provoked by the main body of this work originated 
in the finding that the temperature to which we heat gas 
in AGN feedback is critical. We extend our investigation 
to varying supernova temperature in Section 15.31 Finally, 
we disentangle the acc retio n and feedback e l ement s of the 
iDi Matteo et al.l l|2005h and lBooth fc Schavd l|2009l ) models 
incorporated in the main investigation to confirm that the 
AGN feedback temperature is the main feature driving the 
change in the star formation rate in Section [5.61 



one tenth of the number of particles but with all other galaxy 
parameters held constant (see Table [l}. We begin by probing 
the differences in the star formation rate for simulations run 
at lower resolution. 

The star formation rates in simulations of isolated 
galaxies show a small decrease in star formation due to a 
reduced ability to resolve high-density clumps where stars 
form. However the dependence on resolution is slight, with 
around a ~ 10 per cent reduction, so we do not discuss 
this further. The resolution dependence of simulated merg- 
ers is more complex. Fig. [T^] shows the ratio (low resolu- 
tion/high resolution) of stellar mass formed against time 
for the NFB, SN, KAGN_SN and SAGN^N simulations 
(simulations without supernovae show identical trends). As 
found in isolated galaxies, low-resolution simulations gener- 
ally under-produce stars, leaving more gas. However, they 
show an increase in star formation during starbursts as this 
gas is compressed during the merger. The NFB, KAGN, 
SAGN and DAGN simulations produce 1, 5, 15 and 17 per 
cent less stellar mass over the course of the simulations 
respectiveljQ. The larger reduction in stellar mass for the 
strong feedback AGN models (DAGN and SAGN are very 
similar at lower resolution) is due to the powerful outflows 
interacting with more gas, as a result of the larger smoothing 
lengths in low resolution simulations. 

The black hole mass in low resolution simulations, plot- 



5.1 Numerical resolution 

In order to evaluate the effect of changing the simulation 
mass resolution, we re-run the main suite of simulations with 



* We note that changing the simulation mass resolution does not 
directly modify the time interval or mass loading of feedback 
events for a given accretion rate in the strong feedback method 
because we decrease Afji^at by the same factor by which we de- 
crease the resolution (to A^hcat = !)• 
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Figure 12. Ratio of stellar mass formed during simulations per- 
formed at high and low resolution for NFB, SN, KAGN_SN and 
SAGN_SN models in merging galaxies. The SAGN_SN simulation 
exhibits the strongest resolution dependency due to the difficulty 
in resolving a columnated energetic outflow at low resolution. 
Under-resolving the outflows increases the amount of gas with 
which they interact with and heat. 



ted as a fraction of the high resolution mass is shown in 
Fig. [13] for isolated galaxies (top panel) and mergers (bot- 
tom panel) , for the main AGN simulations. All AGN models 
produce a larger black hole mass in low resolution simula- 
tions. This resolution dependence is due to more poorly re- 
solved accretion, meaning that in low resolution simulations 
a black hole can accrete much more mass before capturing 
a particle and lowering the local density. The effect is even 
more severe in the KAGN model, where the typical heat- 
ing temperature decreases in the low resolution case as the 
kernel mass increases. 

A direct co mparison with the bl ack h ole masses found 
in the work of ISpringel et all l|2005h and lOi Matteo et all 
l|200a) is made possible using our lower resolution simula- 
tions, as their runs used 8 x 10'' particles. In isolated galaxy 
simulations we find a black hole mass of 1.8 x 10® Mq , in rea- 
sonably good agreement with Springel et al. (~ 1 x 10® M©) 
albeit with our different resolution, star formation routine, 
ISM model and disc gas fraction. In merging systems we find 
a final black hole mass somewhat lower (3.2 x 10® Mq) than 
the value obtained i n previous studies (~ 1 — 3 x 10^ Mq; 
ISpringel et al.l l2005l : [Johansson et al.l |2009| ) . Such a devia- 
tion may be due to our choice of irregular m erger geometry 
which is known to reduce black hole mass (|Springel et al.l 

[iooi). 



5.2 Initial gas fraction 

We now investigate how the supernova and AGN feed- 
back processes interact in gas-rich discs with a gas fraction 
/g = 0.9, more appropriate for a high redshift system. (The 
simulations are performed at low resolution for the standard 
AGN models to minimise CPU time consumption.) Fig. 1141 
shows the black hole mass evolution for the low resolution 
fiducial models (left panels) and increased gas fraction (cen- 
tre panels) . This change has the effect of accelerating the co- 
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Figure 13. Total black hole mass ratio in low resolution over 
high resolution simulations in isolated (top) and merging (bot- 
tom) galaxies. For all models, black holes grow significantly more 
in low resolution simulations. 



alescence of the galaxies/black holes. However, we find that 
the black holes in gas-rich systems cease growth at broadly 
similar masses to the fiducial simulations for all AGN mod- 
els. This is perhaps surprising as we have shown that the 
black hole growth is very sensitive to the amount of gas avail- 
able in its surroundings (Section 14. 4p . This is explained by 
the increased level of star formation observed in the gas-rich 
galaxies (Fig. I15|l consuming the gas on the same timescale 
as in the fiducial galaxy and terminating black hole growth. 
Aside from the higher star formation rates seen across all 
times, the feedback models show the same trends in the 
star formation rate as in the case of the models with lower 
gas fraction. Supernovae dominate the star formation sup- 
pression in the isolated galaxies and through mergers the 
AGN model with weak feedback suppresses the star forma- 
tion slightly, but it is shut off completely by the strong feed- 
back models. 

5.3 Initial black hole mass 

The initial black hole mass used in our main simulation 
runs is motivated by the approximate mass scale of pri- 
mordial 'seed' black holes formed th rough the direct col- 
lapse of zero metallicity gas (see e.g. iBellovarv et al.ll201ll 
and references therein) and the desire for compatibility with 
previous studies. This value is poorly constrained and well 
below the mass of a black hole predicted by observational 
relations for this galaxy, Mbh ~ 10^ Mq ( Magorrian et al.l 
1 19981 : iTremaine et al.ll2002l : iHaring fc Rixll2004l ). Here we in- 
vestigate the sensitivity of the black hole mass evolution to 
the initial mass chosen. We perform additional low resolu- 
tion simulations for the KAGN, SAGN and DAGN models 
with an initial mass of A/bh,! ~ 10® Mq. While this mass 
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Figure 14. Black hole mass throughout low resolution, isolated (top) and merging (bottom) galaxy simulations including the three 
primary AGN models. Columns (left to right) show: standard initial black hole mass and disc gas fraction (A/bh,! = 10^ Mq, /g = 0.3); 
increased disc gas fraction (Mbh.i = 10^ Mqi /g = 0-9); S'ld increased initial black hole mass (Mbh i = W^ Mq, /g = 0.3). The final 
black hole mass is relatively unaffected (to within a factor of 2-3) by changing either the disc gas fraction or initial mass. 
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Figure 15. Star formation rate throughout simulations of low 
resolution merging galaxies with an increased disc gas fraction 
/g = 0.9 including the three primary AGN models and super- 
novae (KAGN.SN.g9, SAGN.SN.g9, DAGN.SN.g9). The shape 
and trends are the same as seen with the fiducial gas fraction, 
however the absolute rate is almost an order of magnitude larger. 



Figure 16. Black hole mass evolution for isolated galaxies sim- 
ulated with the fiducial AGN models and no star formation. 
This shows that the black holes have not self-regulated their own 
growth through feedback within 3 Gyr. 



final mass in tlie standard case. Tlie star formation rate is 
similarly independent of initial black hole mass (not shown). 



is substantially larger (10 times the fiducial mass) it is still 
lower than the predicted value and therefore still requires 
significant growth to reach this relation. 

We find that changing the initial mass in this way (Fig. 
I14p does not greatly affect the final mass of the black holes. 
Growth ceases at similar values to the fiducial low resolu- 
tion runs provided that the initial mass does not exceed the 



5.4 No star formation 

We showed previously (Section [3.2p that in isolated galaxies 
the growth of the central black hole slows over time. This 
may be interpreted as evidence that the black hole has self- 
regulated its growth through feedback, however we noted 
that this 'self regulation' occurs on a similar timescale as 
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Figure 17. Projected logarithmic faee-on gas surface density maps for isolated galaxies simulated with different supernova feedback 
parameters at t = 0.6 Gyr (all simulations include Bondi-type kernel-weighted AGN feedback). Models shown have increasing (left to 
right) feedback strength: (KAGN_SN_)N10T0.1, NlTl, N3T1, N1T3. A larger heating temperature leads to an increased porosity of the 
gas disc. 



Table 3. Parameters for supplementary simulations investigat- 
ing supernova feedback. Columns detail the simulation label, su- 
pernova feedback efficiency, number of gas particles heated by 
each star particle and the temperature by which gas is heated 
through a feedback event. Note that KAGN_SNJN[1T1 is our fidu- 
cial choice. 



Label 


eSN 


SN AThcat 


Thoat [K] 


KAGN.SN.N1T1 


0.37 


1 


1 X 


10^ 


KAGN.SN.N10T0.1 


0.37 


10 


1 X 


lO** 


KAGN.SN.N3T1 


1.1 


3 


1 X 


lO'^ 


KAGN.SN.N1T3 


1.1 


1 


3x 


lO'^ 



gas is depleted from the galaxy. As a means of modifying 
the timescale over which gas is consumed without substan- 
tially altering the dynamics or disc structure, we perform 
additional simulations in which we artificially turn off star 
formation. 

Fig. [16] shows black hole mass versus time for models 
with no star formation. All models yield substantially larger 
final black hole masses and the slowing of black hole growth 
towards the end of the simulation is less clear. This shows 
that, at least in artificially isolated galaxies without an ex- 
ternal gas source, black hole growth is halted due to external 
gas consumption through star formation. Although a galaxy 
would normally be fed with cooling gas from a halo or ac- 
cretion along a filament, it is possible that the regulation 
of black hole growth through star formation is important in 
low redshift galaxies for which filamentary structures have 
been disrupted and/or the gas becomes too diffuse to cool 
efficiently onto the galaxy. 



5.5 Supernova feedback parameters 

As a probe of the numerical sensitivity of our supernova 
feedback implementation we also present additional low res- 
olution simulations performed with varying mass loading 
factors and supernova heating temperatures. In the fiducial 
model, each star formed causes one neighbouring gas par- 



ticle to be heated to lO'^K, coupling the supernova energy 
with an efficiency, esN = 0.37. Here we additionally simulate 
this feedback with several different parameter choices (Ta- 
bleO, namely heating 10 particles to 10^ K, thereby keeping 
the total energy output constant, a case where we deposit 
all the energy available (see Section r2.3l for a derivation) into 
3 particles heated to lO'^ K; and a maximal case where we 
heat 1 particle to 3 x 10^ K. 

Fig. [17] shows projected gas density plots for the range 
of supernova parameters. Increasing the strength of the feed- 
back, and particularly the temperature, increases the poros- 
ity in the gas disc and redistributes more gas into the hot 
halo. Note that the simulations were performed with the 
KAGN method in conjunction with supernova feedback as 
this model allows an investigation of the effect on the black 
hole during a merger, without the AGN terminating star 
formation. 

Fig. [18] shows the evolution of the star forma- 
tion rate in isolated (top panel) and merging (bottom 
panel) galaxies. Star formation is barely suppressed in the 
KAGN_SNJv[10T0.1 simulation, whilst the fiducial model 
KAGN_SNJN'1T1 is only marginally less powerful than the 
KAGN_SNJ^3T1 variant, despite the latter having three 
times as much energy available per supernova. Finally, the 
KAGN_SNJN'1T3 model causes by far the most suppression 
of star formation at all times, even delaying the timing of the 
starburst. The delay does not have a dynamical origin as the 
black holes are dynamically bound/merged by i = 1.37 Gyr. 
It is caused by the longer time it takes for the gas to con- 
dense out of the hotter and more diffuse halo created by the 
supernovae. 

In summary, the level of star formation suppression due 
to supernovae is most sensitive to the choice of heating tem- 
perature, a finding in agreement with the results from our 
AGN feedback simulations. 



5.6 AGN feedback parameters 

The findings of our main investigation indicate that the most 
important aspect of an AGN model, in terms of the effect 
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Figure 18. Evolution of the star formation rate with varying su- 
pernova feedback strength for low resolution simulations employ- 
ing kernel-weighted, Bondi-type AGN feedback in isolated (top 
panel) and merging (bottom panel) disc galaxies. The low resolu- 
tion fiducial model (KAGN) is shown for comparison. The N1T3 
model causes by far the most suppression of star formation. 



on the host galaxy, are the details of how the feedback en- 
ergy is deposited. Specifically, the temperature reached by 
heated gas is shown to be critically important. Striking dif- 
ferenc es were found between the Bondi-based AGN mod - 
els of iDi Matteo et al.1 (|2005| ) and iBooth fc Schavd (|2009l ) 
which differ primarily by employing kernel-weighted and 
strong feedback respectively. However, the details of the two 
methods have additional variations; the values of the feed- 
back coupling efficiency (ef — 0.05 versus 0.15) and accretion 
rate boost, a — 100 or oc p^ (see Section [2.4.11 for details) 
both change between the models. In order to disentangle the 
effects of these choices we have performed two extra simula- 
tions, changing each of the parameters one at a time. Table 
[4] summarises the extra simulations performed. 

In Fig. 1191 we compare the star formation rate (top 
panel) and total black hole mass (bottom panel) between 
the fiducial Bondi models and the additional merger simu- 
lations at low resolution. This shows conclusively that the 
heating method is the dominant factor in determining the 
star formation suppression by AGN feedback in these sys- 



Table 4. Additional simulation parameters probing the AGN 
feedback method. Columns detail the simulation label, Bondi ac- 
cretion rate boost factor, number of feedback-heated particles, 
feedback distribution type and AGN feedback efficiency. 



Label Ace. 
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Type 


tf 
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Figure 19. Effect of varying AGN feedback parameters on the 
star formation rate (top panel) and total black hole mass (bottom 
panel) in low resolution merger simulations. The star formation 
rate through a merger is dominated by the choice of feedback 
model, while the final black hole mass also shows a strong depen- 
dence on the efficiency parameter Cf. 



tems. Changing the accretion rate density dependence does 
have a small effect on the star formation at early times, but 
the dominant trend of declining and terminating star forma- 
tion still occurs at the same time. Modifying the available 
energy in feedback by factor of three through the et param- 
eter has no discernible impact up on the star formatio n rate 
(a result similar to that found bv lDebuhr et al.ll201ll ). 

The black hole mass evolution is also sensitive to 
changes made to the heating efficiency parameters, showing 
a trend for reducing mass with increasing feedback strength 
as would be expected. These findings show that the black 
hole feedback efficiency can be used to tune the black hole 
mass without affecting the host galaxy, as lon g as the feed- 
back m ethod is constant, in agreement with iBooth fc Schavd 
(|2009l ). 
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6 CONCLUSIONS 

AGN feedback models are becoming a standard component 
in both large scale A'^body hydrodynamical simulations and 
semi-analytic modelling, but the full implications of incor- 
porating even the current simplistic models are not yet fully 
understood. One particular uncertainty is how the combina- 
tion of AGN and supernova feedback processes change the 
star formation rate in a galaxy. We performed and compared 
a series of identical high resolution disc galaxy simulations 
with a range of sub-grid model combinations including a 
variety of AGN models and supernova feedback strengths. 
Changing only the feedback models allowed us to perform 
a clean comparison between competing parametrisations, as 
well as probing the details of the processes involved and 
how the models interact. Beginning with an isolated disc 
galaxy we investigated how feedback processes change the 
star formation rate and where appropriate the black hole 
growth. Building on this understanding we simulated the 
galaxies undergoing major mergers and analysed the effects 
of AGN and supernova feedback on the galaxy through an 
active phase in its evolution. Our main conclusions may be 
summarised as follows: 

• In isolated Milky- Way-mass model disc galaxies super- 
nova feedback dominates the star formation suppression, 
whereas AGN feedback has no impact for any of the tested 
models. This is due to the disparity in scales upon which 
the processes act; AGN activity is (eponymously) confined 
to the central region whereas supernovae occur throughout 
the galactic disc. 

• AGN feedback plays a much larger role in mergers and 
may completely eliminate the starburst when the heating is 
strong. The most important factor determining the strength 
of AGN feedback on the host galaxy is the temperature to 
which gas is heated and correspondingly the feedback mass 
loading. A higher temperature allows (fewer) gas particles 
to escape from high density regions without suffering from 
significant radiative losses. 

• Supernova and AGN feedback are largely independent, 
with only a slightly stronger suppression of the star for- 
mation rate when acting together than if they acted in- 
depen dently. This finding is in contradiction to the find- 
ings of iDalla Vecchia fc Schavd (|2012i ) an d supports the ap- 
proach taken in semi-analytic model ling, (IBower et al.luOOq : 
ICroton et al.l 120061 : IGuo et al.l 120111 ). However, our simula- 
tions are a restricted set of conditions and ignore complex 
effects such as external gas accretion. 

• Numerical resolution strongly affects the black hole 
mass due to under-resolved feedback and accretion. The 
problem is more severe for the kernel-weighted model where 
lower resolution simulations result in lower heating temper- 
atures due to the increase in mass within the kernel. How- 
ever, the stellar mass shows the most resolution dependence 
in strong feedback models. 

• The black hole mass is sensitive to the choice of both 
the feedback and accretion models, however the parameters 
may be tuned to obtain desired masses with little effect on 
the galaxy. The mass is strongly affected by gas consumption 
through star formation in these model systems which neglect 
an external gas source such as a hot halo or minor mergers. 
Gas fraction and initial black hole mass do not affect the 
final black hole mass by a large factor. 



• The level of star formation suppression due to super- 
novae is (as for AGN feedback) strongly linked to the tem- 
perature of heated gas with higher temperatures causing a 
larger reduction. The number of heated particles (mass load- 
ing) is of secondary importance. However, even the maxi- 
mally strong supernova feedback model was unable to pre- 
vent a starburst occurring in the galaxy merger when AGN 
feedback was weak. 

While we have gone some way towards disentangling 
the role of feedback processes across multiple models for 
the same systems, further work is required to make progress 
with many aspects of our understanding of AGN feedback in 
galaxies. Due to our need to perform multiple high resolution 
simulations, we were only able to study 2 identical systems. 
Future work, taking advantage of advances in computational 
power, may i nvestigate these pro cesses in cosmological zoom 
simulations (|Tormen et al.lll997l ) of multiple galaxy masses 
and environments. A possible shortcoming of this work is 
that the IMF employed (Salpeter) is currently in question, 
with alternatives e.g. Chabrier and Kroupa IMFs arguably 
giving better fits to data. An alternative IMF alters the star 
formation history of a galaxy as well as the energy available 
from supernovae. However the trends outlined in this paper 
are unlikely to change substantially with an alternative IMF 
for a given supernova model. 

Despite ongo i ng work (e-K - | Hobbs et al.l |2012| : 
iDebuhr et all I2OIII : iFanidakis etall I2OIII : iPlaneUej |2012| ). 
the current generation of AGN models such as those 
incorporated in this study are unable to perform well over 
the full mass range in cosmological objects, from dwarf 
galaxies to clusters, as a fully developed model should. It 
is therefore important that the development of a flexible, 
resolution independent AGN feedback model continues. 

While this paper was in advanced dra ft form, two recent 
papers appeared by IWurster fc Thacken (|2013al lbl). These 
papers share some commonality with this work as they 
also investigate several AGN models in idealised isolated 
and merging disc galaxies, including the 1 Power et al.i (|201ir ) 
disc accretion model. The papers are complementary to our 
own as they choose to focus primarily on the details of the 
AGN models themselves, rather than their interaction with 
the host galaxies. A stronger star formation suppression is 
found for strong AGN feedback than kernel-weighted feed- 
back (in their corresponding models), in agreement with our 
own findings. 
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APPENDIX A: CREATING EQUILIBRIUM 
DISC GALAXIES 

In creating our model galaxies we largely follow 
ISpringel et al.l ((20051); however our approach differs in a few 
key areas such as our treatment of adiabatic contraction of 
the DM halo (see Section IA1.3|) and the ISM model used 
(Section 12. 2p . The model galaxies created here will consist 
of an encompassing DM halo surrounding a disc made up of 
stars and gas particles with an optional stellar bulge. 



Al Collisionless particle positions 

When defining the positions of the collisionless components 
we begin with a so-called 'glass', a volume of random and 
homogeneously distributed particles. We then modify the 
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glass by adding additional constraints on the distribution 
of the particles to fit desired profiles. In this way any unde- 
fined aspect of the particle distribution is inherently random 
and we do not rely solely on pseudo-random number gener- 
ators. Should the glass in use contain insufficient particles 
we copy and tile the existing glass until we obtain a cube 
large enough. 



Al. 



Dark matter halo & the stellar bulge 



For the DM comp onent we a s sume that it follows the Hern- 
quist mass profile iHernauisj l| 199(3 ) 



Ph(r) 



M, 



27r r(r + a)^ ' 



(Al) 



where Mh is the total halo mass (A/h — fhM where M is the 
total galaxy mass and /h is the halo mass fraction) and a 
is a constant scale length which affects the central shape of 
the profile. The corresponding enclosed mass profile is then 



Mh(<r) = A4- 



(r-l-a)^ 

We choose this profile as it replicates t he inner slope of the 
cosmologically motivated NFW profile (JNavarro et al.lll996l ) 
whilst having a mass profile which converges at large radii. 
We utilise the following conversion between a and the NFW 
concentration parameter, c, under the assumption that the 
two profiles contain the same mass within r2oo (the radius at 
which the density of the system drops below 200 x the crit- 
ical density required to close the universe at that redshift) 



a^Ts V2[ln(l+c)-c/(l + c)], (A3) 

where Vs = r2oo/c is the scale length of the NFW halo. We 
then distribute the DM particles simply by radially sliding 
the particles to the position found by interpolating Equation 
[A2l at M{< r) = {N + l)mDM, where iV is the number 
of particles already placed and nioM is the mass of a DM 
particle. In this way the random angular position of the 
particle is maintained. 

For the stellar bulge we also follow a He rnquist profile 
with the modification that a different flarger. lNavarro et al.l 
Il99a ) scale length b is employed instead of a and the particle 



mass IS now nistar 



Mi,{<r) = A-h- 



(A4) 



' (r + by • 

Both a and b are free parameters of our model. 

A1.2 Stellar disc 

For the stellar disc we define that its surface density is 
described by an axisymmetric exponential disc with scale 
length h 



E.(i?) 



A-R/h) 



27r/i2 



(A5) 



where Af*(= (1 — fg)Md) is the mass of the stellar disc and 
/g is the disc gas mass fraction (Md = M« -|- Mg). From this 
we deduce the enclosed mass profile and perform the same 
procedure as followed for the halo when setting the particle 
radii, sliding in R 



M.«i?) = M.^^±:^e(-«/''). 



The scale length, h, is set for both gaseous and stellar com- 
ponents by relating it to the angular momentum of the disc. 
We calculate the angular momentum by assuming that the 
baryonic matter which now makes up the disc was initially 
distributed identically to the surrounding DM halo forming 
a 'total' halo and then applying conservation of angular mo- 
mentum. In this case the angular momentum of the disc is 
equal to a fraction of the 'total' halo angular momentum, 
i.e. 

Jd = id J, (A7) 

where we take the approach of ISpringel fc White! (| 19991 ) in 
assuming there is no angular momentum transport between 
the different galaxy components, under this assumption the 
angular momentum fraction, (jd, is equal to the mass frac- 
tion of the disc. The spin of a halo is often described by the 
dimensionless spin parameter, A 

where E is the total energy of the halo. For an NFW halo 
the kinetic energy is found to be, under the assumption that 
all particles move on cir cular orbits at the circular velocity 
jSpringel fc Whitdl 19991 ) 

f2 



where 



/c = 



-Bkin 



c 1- 



2r2oo 



fc, 



(l+cV 



21n(l + c) 
1 + c 



2[ln(l + c)-^ 



(A9) 



(AlO) 



Although this relationship is strictly true for NFW profiles 
we assume it also holds for the Hernquist profile due to the 
similarity of the profiles in the central regions. Then apply- 
ing the Virial theorem so E = — -Ekin and by substituting 
into Equation IA8I we obtain an expression for the halo an- 
gular momentum 



J = XG^'^M. 



1/2,^3/2^1/2 



(All) 



We may then find a h value for which the numerically cal- 
culated Jd value agrees with the predicted value (in prac- 
tice DM halo contraction must be simultaneously consid- 
ered, see Section rA1.3[) . The angular momentum of any disc 
with known rotation curve and surface density may be found 
with (JMo et al.lll9"9i ) 



Jo 



Jd = Mi VciR)R X 2nRT,{R) dR, (A12) 



where 14 is the circular velocity. For our model we substitute 
in Equation IA5I to find 

Jd = Ma r Vc{R) {^] expf^^^dR (A13) 

where Afd = is the total disc mass, and the circular velo city 
in our system is given by (e.g. iBinnev fc Tremaindl 19871 pg. 

77), 

,2,r>^ G[M^{< R) + M^{< R)] 



V,\R) = 



(A6) 



R 

+ '^^y^[h{y)K,{y) - h{y)K,{y)], (A14) 
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where y = R/{2h) and /„, Kn axe Bessel functions. 

For the vertical distribution we employ the model of an 
isothermal sheet with constant scale height 20 ■ From this 
we once again find the enclosed mass profile which can be 
rearranged to give the vertical positions 

S.(.) = ^sech^ (^\ , (A15) 



2zo 



M^{< z) = M^tanh 



2zo 



(A16) 



Here we follow ISpringel et al.1 l|2005l ) in setting zq = 0.2/i. 
For completeness, the 3D density and mass distribution of 
the stellar disc are described by 



(R z)- — M^sech^ ( — \ e'-^/'^' 



(A17) 



(^)(.-e--(^ + .)) 



M*(< i?, < 2) = — -tanh 

(A18) 

We note here a numerical consideration that the R, z 
coordinates of a particle are independent, therefore any ar- 
tificial correlation must be carefully avoided. 

A 1.3 Halo contraction 



As an optional feature we include the ability to approximate 
the DM halo contraction resulting from the formation of a 
•alact ic disc. Following the approach of iBlumenthal et al.l 
19851) (also IMo et al.l Il998l ) we assume that the matter 
which makes up the galaxy, both dark and otherwise, was 
initially distributed following a Hernquist profile (as em- 
ployed in our calculations of the disc angular momentum) 
but as the baryonic matter began to cool and collapse the 
DM distribution responded by contracting adiabatically. We 
also assume that the halo remains spherically symmetric un- 
der the contraction, under these assumptions the angular 
momentum of individual DM particles is conserved. Consid- 
ering a single particle at initial radius, ri, moving to a new 
position with radius, rt, we may write 



GMt{< n) = GA/i(< n). 



(A19) 



where Mi is the total initial mass enclosed within ri as given 
by the Hernquist profile and M{ is the total final mass en- 
closed within Tf. The final mass Mt is equal to the sum of 
the (now contracted) DM component of the mass initially 
within ri {= /hAfi(< ri)) and the baryonic matter within rt 

M{{< n) = Md(< rt) + Afb(< rf) + /hAfi(< n), (A20) 



where we analytically know Afb(< r) from Eg nation I A4I and 
Md{< r) is found similarly to Equation IA18I whilst setting 
R = z — r 

Afd 



A^d(< r) 



■ tanh , , 

h \ 2zo ) 



fe ('-"'"""■ + '") 



(A21) 



By iterating we then find a value of rf which satisfies Equa- 
tion IA19I The problem is further complicated because the 
determinations of the disc scale length, /i, and the halo con- 
traction are interdependent; we must therefore iterate over 
the two calculations from a first guess until the solutions 
converge. See Fig. lAll for an illustration of the effects of 
contraction on rotation curves. 
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Figure Al. Rotation curves as derived with (top panel) and 
without (bottom panel) adiabatic contraction of the DM halo. 
Also shown as a vertical dotted line is 2.8e where e is the Plummer 
equivalent softening length. Note that the two cases have differing 
disc scale lengths, h = 3.45 kpc and 2.78 kpc respectively, this is 
due to the dependence of the scale length on the DM halo circular 
velocity and therefore its contraction. 



A 1.4 Gaseous disc 

The radial dependence of the gas distribution is set as for 
the stellar component in Section TA 1.2 1 and follows 



E,(i?) = 



M, 



27rft2 



^g{-fl/h)^ 



(A22) 



To ensure the gas is initially placed in a stable structure 
we assume hydrostatic equilibrium 



dP a$ 

= -Pg- 



d. 



dz 



(A23) 



where $ is the total gravitational potential. Substituting in 
our EoS (Equation m we obtain 






(A24) 
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Rearranging we see 

^^ Az A-y Az ' 

separating variables then gives us the integrals 



„7-2 



'^P- = -T^ 



d$, 



(A25) 



(A26) 



where po,"l'o are the respective mid-plane values. Performing 
the integrations and rearranging gives 



Pgiz,R) = 



Po 



'\R)- 



1 



A^ 



mz,R)-MR)] 



(A27) 

The mid-plane potential <I>o is known numerically; however 
we must determine the mid-plane density, po, by iteration 
demanding that the 3D density agrees with the surface den- 
sity i.e. that it fulfills the constraint, 



Sg(i?)= / pg(i?,2,$)dz. 



(A28) 



At a given radius we may then obtain a profile for Mg{< z) 
by numerically integrating Equation I A27I we then interpo- 
late this profile when vertically positioning the gas particles. 
The known density distribution is then used with Equation 
[T]and the ideal gas law to find the temperature of each par- 
ticle. 

At this stage we discuss an important numerical consid- 
eration in this method, many of the preceding derivations 
have required the use of spatial derivatives and we have be- 
gun to introduce the use of the potential, $, for both of 
these purposes we create a fine logarithmic mesh in R, z. It 
is at these mesh points which we calculate the gravitational 
potential due to the mass distribution using a parallelised 
'Tree' method and define the density profiles. The mesh is 
also employed extensively when calculating the velocity dis- 
tributions. 



A2 Particle velocities 

A2.1 Collisionless particle velocities 

When considering the velocities of t he particles represent- 
ing collisionless components we follow iHernguistl (|l993l ) (see 
also iBinnev fc Tremaind 119871 ) and calculate the distribu- 
tion by considering the collisionless Boltzmann equation 
(CBE). In deriving our distribution we assume that it only 
depends on the energy, E, and the ^-component of the an- 
gular momentum Lz as well as assuming that the velocities 
are isotropic. In such a case the mixed second order mo- 
ments vanish as do the radial and vertical first order mo- 
ments ((iiRUz) = {vy{.v4,) = {vzv^) = 0, {vn) = ("z) = 0). 
The non- vanishing second order moments are then obtained 
from the Jeans equations as 



(«z) = (^r) = - / p{z,R)-—Az, 
P Jz oz 



{■"I) 



K) + 



R d{p{vi)) 
p dR 



(A29) 



(A30) 



where p is the density of the currently considered mass com- 
ponent, $ is the gravitational potential from all components 
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Figure A2. Evolution of galaxy rotation curves for an adiabat- 
ically contracted isolated model galaxy with 7 X 10® particles. 
Rotation curves are shown at 4 = 0.0 and 3.0 Gyr in a grav- 
ity and hydrodynamics only simulation for (top to bottom) dark 
matter, stellar disc, gas disc and stellar bulge. Also shown as a 
vertical dotted line is the estimated two-body relaxation radius 
at 1.0 Gyr. A small level of contraction in the central region is 
seen, however for the most part the galaxies are extremely stable. 



and Vc is the circular velocity 



R 



dR' 



(A31) 



We then require a suitable function to approximate the ve- 
locity distribution at any point, for this purpose we choose a 
triaxial Gaussian distribution with its axes coincident with 
the axisymmetric coordinate system 



F{vx,R,z) = 



: exp 



{vx - pY 
2a\ 



(A32) 



where X(= R,z,(f>), crx(= '^x{Rj^)) is the velocity vari- 
ance along a general axis, X, and p{— p{R,z) = {vx}) is 
the mean streaming velocity. The velocity dispersions and 
variances are related by 

2 / 2 1 

ai = (vl) (A33) 

f^l = {i^l) - M^ 

where we have used (ur) = (vz) = 0. 

For the stellar bulge w e set t he streaming velocity, p, 
to zero as in ISpringel et al.l (120051 1 however this is not true 
for the DM halo. For the halo we set the mean rotational 
velocity at any point to be a fixed fraction of the circular 
velocity (with a fiducial value of /i = O.Oluc), note however 
that this is of small importance when compared to the radial 
motions of the DM particles and it is the velocity dispersion 

which contributes the majority of support to the halo. 

Once more we take the approach of ISpringel et al.l 

(|2005|) and continue to employ the Gaussian approximation 
when describing the stellar disc with the assumption that 



("r) = Mivz), 



(A34) 



where it is assumed that /r = 1, this is justified as being 
well supported by observational evidence. In the case of the 
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disc, unlike the halo and bulge, it is the rotational motions 
of the particles which is most important for providing sup- 
port to the structure. When considering this mean azimuthal 
streaming of the disc we invoke the epicyclic approximation 

2 



_2 2 K_ 



^4- = ^^l^' (A35) 



where the epicyclic frequency, k is defined as 
2 _ 3_a£ 9^ 



+ 77^' (A36) 



and 



We then define 



" = ^i- (^^^) 



to obtain the simplified relation, a^ — a'^/rf. Equations 
IA30I and IA33I then allow us to write the mean streaming 
velocity 



("0> = a/W>-^- (A39) 

Numerically, the velocity distributions are defined for 
the fixed spatial mesh described briefly in Section [A 1 .41 once 
the matter distributions have been finalised. This is so that 
the full gravitational potential is then known in regularly 
distributed positions and then may be integrated or differ- 
entiated as needed. When finding the velocity of any given 
particle, bilinear interpolation is used to find the approxi- 
mate a (and where appropriate y,) at the particle's position, 
the Gaussian distribution is then sampled using the so-called 
accept-reject method. 



A2.2 Gaseous disc velocities 

The velocity structure of the gaseous disc is much simpler 
because the majority of its support is provided by its pres- 
sure, it is therefore assumed to undergo solid body rotation 
i.e. its velocity depends only on R and vb. = vn = 0. As- 
suming balance between the gravitational pull and the cen- 
trifugal and pressure support on the disc we have 



9$ _ «^,gas 1 dP 



(A40) 



dR R pgdR' 

where P and pg are the gas pressure and densities respec- 
tively. We then find the particle velocities to be 



,„/a$ 1 dP\ 

The model galaxy is now fully specified and from the 
method and assumptions we have utilised we expect it to 
show very little evolution over multiple dynamical times, as 
shown in Fig. 
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